Untitled

Introduction

The tbl_summary() %>% add_p() function from the gtsummary package is updating the default test for most continuous variables. Previously, continuous variables being compared between two groups would be compared with chisq.test(correct = TRUE), which would use a continuity correction for 2 by 2 table comparisons. However, the simulation study below demonstrates the superior properties of chisq.test(correct = FALSE), and we have updated the default to use this test instead.

Simulating results

We will simulate many studies assuming various prevalence rates of two binary variables and at various sample sizes. The power and size of the chi-squared tests will be reported.

Attach needed libraries.

library(tidyverse)
library(gtsummary)

Set the number of simulated studies that will be run for each combination of rates and sample size.

simn <- 100

Simulate the trial results.

sim_results <-
  list(
    simn = seq_len(simn),                              # number of simulations per scenario
    n = c(20, 50, 100, 500),                           # sample sizes
    p = c(0.05, 0.10, 0.20, 0.50),                     # rate of success in group 1
    delta = c(0.0, 0.05, 0.10, 0.20, 0.30, 0.40, 0.45) # p + delta = rate of success in group 2
  ) %>%
  cross_df() %>%
  # deleting duplicated scenarios (due to symmetry of test)
  filter(p + delta <= 0.5) %>%
  mutate(
    # simulating  data set with specified success rates, and calculating p-values
    df_result = pmap(
      list(n, p, delta),
      function(n, p, delta) {
        data = 
          bind_rows(
            tibble(value = runif(n/2) < p, group = 1),
            tibble(value = runif(n/2) < p + delta, group = 2),
          )
        
        tribble(
          ~test, ~reject, 
          "With Correction", with(data, table(value, group)) %>% 
            {tryCatch(chisq.test(., correct = TRUE), error = function(e) NA)$p.value} %>% 
            {. < 0.05},
          "Without Correction", with(data, table(value, group)) %>% 
            {tryCatch(chisq.test(., correct = FALSE), error = function(e) NA)$p.value} %>% 
            {. < 0.05}
        )
      }
    )
  ) %>%
  unnest(df_result) %>% 
  # deleting sets of results where at least one of the p-values are missing
  group_by(simn, n, p, delta) %>% 
  filter(sum(is.na(reject)) == 0) %>%
  ungroup()

Summarize simulation results.

# calculating power (or level) of test across simulations
results <- 
  sim_results %>%
  group_by(n, p, delta, test) %>%
  summarise(power = mean(reject),
            simn = n(),
            .groups = "drop") 

Power of test

Across sample sizes and various prevelances of the binary variables, teh test without continuity correction exhibits higher power.

results %>%
  filter(p != 0.5) %>%
  mutate(n = str_glue("N = {n}") %>% factor(levels = unique(.)), 
         p = str_glue("{p * 100}%") %>% factor(levels = unique(.))) %>%
  ggplot(aes(x = delta, y = power, color = test)) +
  geom_line() +
  geom_point() +
  geom_hline(yintercept = 0.05) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
  facet_grid(n ~ p, scale = "free") +
  labs(color = "",
       y = "Power",
       x = "Delta") +
  theme(legend.position="bottom")

Level of test

In the simulation study, we rejected the null hypothesis when the p-value was less than 0.05. The true rate of rejection should be 5%, and the test without continuity correction is closer to this threshold.

results %>%
  filter(delta == 0) %>%
  mutate(n = str_glue("N = {n}") %>% factor(levels = unique(.))) %>%
  ggplot(aes(x = p, y = power, color = test)) +
  geom_line() +
  geom_point() +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
  facet_wrap(~n) +
  geom_hline(yintercept = 0.05) +
  labs(color = "",
       y = "Prob. of Rejecting H0",
       x = "Event Rate") +
  theme(legend.position="bottom")

Restore Previous Behavior

If a user wishes to continue using the continuity correction, there are two options.

1. Specify the variables you’d like to compare with the chi-squared test without continuity correction. For example, this line of code specifies all categorical variables are compared with "chisq.test" that uses the continuity correction: trial %>% tbl_summary(by = trt) %>% add_p(all_categorical() ~ "chisq.test")

2. If you want to make the global change for the entire script or R session, use the gtsummary themes. Run the following line anywhere in your R script, which will make the chi-sequared test with continuity correction the default: set_gtsummary_theme(list("add_p.tbl_summary-attr:test.continuous_by2" = "chisq.test"))